A story of viral co-infection, co-transmission and co-feeding in ticks: how to compute an invasion reproduction number

With a single circulating vector-borne virus, the basic reproduction number incorporates contributions from tick-to-tick (co-feeding), tick-to-host and host-to-tick transmission routes. With two different circulating vector-borne viral strains, resident and invasive, and under the assumption that co-feeding is the only transmission route in a tick population, the invasion reproduction number depends on whether the model system of ordinary differential equations possesses the property of neutrality. We show that a simple model, with two populations of ticks infected with one strain, resident or invasive, and one population of co-infected ticks, does not have Alizon’s neutrality property. We present model alternatives that are capable of representing the invasion potential of a novel strain by including populations of ticks dually infected with the same strain. The invasion reproduction number is analysed with the next-generation method and via numerical simulations.


Introduction
Co-infection of a single host by at least two distinct viruses provides an opportunity for viruses to exchange genetic information through genomic reassortment or recombination [1,2].In fact, entirely novel pathogenic viruses have emerged from reassortment events of less pathogenic parents in nature [3][4][5].Co-infection can be thought of as the rate-limiting step in the sudden emergence of genetically distant variants of existing human pathogens such as influenza, SARS-CoV-2, and Crimean Congo Hemorrhagic Fever Virus (CCHFV).Therefore, understanding the dynamics of co-infection in common host species, e.g., arthropods (ticks or mosquitoes), is essential to study the emergence and re-emergence of both new and old human pathogens.
Genomic reassortment is possible in viruses with segmented genomes, such as the Bunyaviruses, which themselves include lethal pathogens of relevance to public health and of pandemic potential, e.g., Lassa fever, Rift Valley fever and CCHF viruses [6].Fig. 1 illustrates the dynamics of reassortment at the cellular level for Bunyaviruses, or more generally for a tri-segmented virus.CCHFV is a tick-borne Bunyavirus, with the potential to reassort, and an increasing geographical range due to the changing climate [7,8].Understanding how adaptable to different hosts this potentially fatal human pathogen is, what role co-infection (as a first step to genomic reassortment) will play in the generation of potential new viral strains, and how those variants will spread among already infected ticks, is a challenge for theoretical biology.
Due to the ability of ticks to carry multiple viruses or viral strains, epidemiologists have considered co-infection in tick-borne diseases [9][10][11][12][13][14][15][16][17], superinfection [9], and co-transmission [18][19][20].Specifically, epidemiologists are interested in quantifying (or understanding) the invasion potential of a novel virus or strain [21,22], given the endemic ability of a resident one [23].In the same way as R 0 provides conditions for successful establishment of a single virus in a susceptible population, the invasion reproduction number, R I , lends threshold conditions for successful invasion of a second virus when the population is endemic with the first one [24,25].For instance, Gao et al. developed a Susceptible-Infected-Susceptible (SIS) model for tick and host populations, and conducted a systematic analysis of invasion by the second virus [18].More recently, Bushman and Antia have developed a general framework of the interaction between viral strains at the within-host level [26].Pfab  Co-infected cells have the potential to generate new viral progeny, different from that of the parental strains: V C ̸ = V A and V C ̸ = V B .B) A co-infected cell can lead to reassortment events, and produce up to 2 3 different reassortants.
endrick [27] for two pathogens [28].Rovenolt and Tate [29] have developed a model of co-infection to study how within-host interactions between parasites can alter host competition in an epidemic setting.Thao Le et al. [30] have studied a two-strain SIS model with co-infection (or co-colonisation) which incorporates variation in transmissibility, duration of carriage, pair-wise susceptibility to co-infection, co-infection duration, and transmission priority effects.Finally, Saad-Roy et al. [31] have considered super-infection and its role during the the first stage of an infection on the evolutionary dynamics of the degree to which the host is asymptomatic.
In the case of plant pathogens, recent experimental studies have shown the complex nature of vector-virus-plant interactions and their role in the transmission and replication of viruses as single and co-infections [32].Allen et al. modeled the transmission dynamics of viruses between vectors and plants, under the assumption that co-infection could only take place in plants [24].Chapwanya et al. developed a general deterministic epidemic model of cropvector-borne disease for synergistic co-infection [33].Miller et al. have shown that mathematical models on the kinetics of co-infection of plant cells with two strains could not adequately describe the data [34].
Current mathematical models of co-infection need to be put in perspective, as previously discussed by Lipsitch et al. and Alizon [35,36].Alizon compared different models of co-infection and raised an issue of nonneutrality [36].He noticed that certain models of co-infection lead to an invasion reproduction number which does not tend to one, in the limit when the invasive and the resident pathogens are the same.To solve this problem, Alizon proposed an alternative model structure, which includes a population of dually infected individuals with the resident pathogen, to achieve the desired neutral invasion reproduction number [36,37].
In this paper, we first present a mathematical model of a single vectorborne virus to understand the role that different transmission routes play in the dynamics of the infected populations.Then, we study the dynamics of two different viruses, or viral strains, in a tick population making use of a classic co-infection model.After performing an invasion analysis, we explain the issue of non-neutrality of the invasion reproduction number, and propose five neutral alternatives.We conclude the paper with a summary of our alternative proposals, their applicability and limitations.

Mathematical model of a single viral strain in a population of ticks and their vertebrate hosts
We consider a tick population feeding on a population of vertebrate hosts, where both populations are susceptible to infection with virus V 1 .The host and tick populations are divided into susceptible and infected subsets.In what follows the number of susceptible hosts (ticks) is denoted n 0 (m 0 ), and the number of infected hosts (ticks) is denoted n 1 (m 1 ), respectively.
The mathematical model considers immigration, death, viral transmission and recovery events in the populations; namely, susceptible hosts and ticks immigrate into the population with rate Φ H and Φ T , respectively.Susceptible and infected hosts can die with per capita rates µ 0 and µ 1 , respectively, whereas susceptible and infected ticks are characterised by the per capita death rates ν 0 and ν 1 , respectively.We assume an infected host can infect a susceptible tick with rate γ 1 , and an infected tick can infect a susceptible host with rate β 1 .Both of these transmission events involve a tick feeding on a vertebrate host, and are referred to as systemic transmission events [38].The virus can additionally be transmitted from an infected tick to a susceptible one via co-feeding [39,40].This occurs when ticks feed on a host in clusters, and close to each other; that is, on the same host and at the same time.In this instance the virus is transmitted by infected tick saliva, with this route of transmission referred to as non-systemic [38].We denote by α 1 the rate at which an infected tick can infect a susceptible one via co-feeding.We assume that transmission events follow mass action kinetics.For example, in the case of co-feeding, and with m 0 and m 1 the number of susceptible and infected ticks, respectively, the rate of infection for the susceptible population is α 1 m 0 m 1 .Finally, once a tick contracts the virus, it remains infected for life [41].On the other hand, vertebrate hosts are characterised by shortlasting viremia [41,42].We, thus, assume that hosts clear the virus with rate φ 1 [43,44].The above set of events are brought together in the following system of ordinary differential equations (ODEs), which describe the dynamics of susceptible and infected hosts and ticks: We note that this system of ODEs (1) has a virus-free equilibrium (VFE), (n ⋆ 0 , 0, m ⋆ 0 , 0), given by

Basic reproduction number
The basic reproduction number, R 0 , measures the mean number of new infections produced by an infected individual (during its lifetime) in a population at the virus-free equilibrium; that is, when the population is completely susceptible [45].R 0 (for the mathematical model (1)) can be calculated making use of the next-generation matrix method [46] as follows.The sub-system of differential equations for (n 1 , m 1 ) is linearised at the VFE, and its Jacobian, J, is then written as J ≡ T +V , with T the 2×2 matrix of transmission events which accounts for new infections in the susceptible population, and V ≡ J − T , the 2 × 2 matrix tracking the changes in the state of the infected populations [46].The next-generation matrix is defined as K ≡ T (−V ) −1 , and the basic reproduction number, R 0 , is given by the largest eigenvalue of K [46].For our system we have with so that which in turn implies If R 0 < 1, the VFE is stable, and if R 0 > 1, it is unstable.The basic reproduction number can be rewritten as where we have introduced the following type reproduction numbers [47] , which represent the contribution of each route of transmission, tick-to-tick, tick-to-host and host-to-tick, respectively, to the total number of new infections (of ticks and hosts) in the susceptible population.R T T , R T H , and R HT correspond to the entries of the next-generation matrix K (see Eq. ( 5)).R HH = 0, since the virus cannot be directly transmitted from an infected host to a susceptible one.The expression of the basic reproduction number for a single virus (see Eq. ( 7)) clearly shows that co-feeding represents a singular route of transmission, compared to systemic routes.For example, β 1 (or γ 1 ) can be very large, but if γ 1 (or β 1 ) is negligible, the contribution to R 0 of viral systemic transmission will be negligible.Therefore, co-feeding events (as characterised by the parameter α 1 ), might maintain an epidemic if R T T > 1.On the other hand, systemic transmission requires both tick-tohost and host-to-tick transmission routes to be non-vanishing, so that there is a chance for R 0 > 1, since as soon as either β 1 or γ 1 are equal to zero, R 0 = 0 in the absence of co-feeding.We conclude this section mentioning a novel network approach (developed in Ref. [48]), to compute the parameters α 1 , β 1 , and γ 1 from first principles.It is reassuring to note that this approach leads to a next-generation matrix with the same structure as K in Eq. (5).

Parameter values
We make use of recent literature to obtain parameter values for the ODE system (1).Table 1 contains a description of each model parameter, together with its plausible ranges and units.Since infection with CCHFV [42,49] or Borrelia [48] is asymptomatic in ticks and vertebrate hosts (but unfortunately not in humans), we assume it does not affect their death rates; that is, µ 0 = µ 1 ≡ µ and ν 0 = ν 1 ≡ ν [18].Given the narrow ranges in Table 1 for Φ H , Φ T , φ 1 and µ, we fix these parameters as follows: Φ H = 1 host per day, Φ T = 2 ticks per day, φ 1 = 1/6 per day, and µ = 10 −3 per day.We derive plausible ranges for the other model parameters making use of Ref. [48], as illustrated in detail in Appendix A.

Visualization of the basic reproduction number
We illustrate the dependence of R 0 on the transmission parameters α 1 , β 1 and γ 1 , Fig. 2, making use of (6), and the parameter values from Section 2.2.Lighter colours correspond to greater values of R 0 (scale on right).Black lines represent a basic reproduction number equal to its critical value of 1.We set Φ H , Φ T , φ 1 , and µ to the values specified in Section 2.2, and set ν = 10 −2 per day (see Appendix A).We consider a different value of α 1 in each panel: on the left, α 1 = 10 −6 , in the middle α 1 = 2 × 10 −5 , and on the right α 1 = 10 −4 (units as provided in Table 1).The corresponding values of R T T are R T T = 1.2 × 10 −2 , R T T = 0.4, and R T T = 2. Along the x-axis and y-axis we vary γ 1 and β 1 , respectively, from 0 to their maximum value listed in Table 1.We note that the area under the curve R 0 = 1 becomes smaller as α 1 increases (from left to right), until it becomes zero when cofeeding transmission contributes to make R 0 greater than one on its own.As one would expect, smaller values of the transmission parameters α 1 , β 1 and Parameter Event Range Units Reference Transmission of one copy of Death rate of T c 10 −2 1/day [48] Table 1: Model parameters introduced in (1) (top half), and ( 8), ( 12), (13), and ( 15) (bottom half).
γ 1 correspond to lower values of R 0 (purple regions on the left and middle panels).Finally, we also note the symmetric role of β 1 and γ 1 in R 0 , as shown in (6).

Two viral strains: tick population and co-feeding transmission
In the previous section, we have shown that co-feeding can sustain an infection among ticks without systemic transmission.The remainder of the paper will focus on the co-feeding route of transmission.We now move to the more complex case where multiple viral strains co-exist, introducing the notions of co-infection and co-transmission.The population of ticks can be infected by two different circulating viral strains, V 1 and V 2 .V 1 is considered to be the resident strain and V 2 the invasive one (e.g., one that emerges once the tick population is endemic with V 1 ).The population of ticks can be classified by its infection status in four different compartments, as susceptible and infected ticks with the resident strain, m 0 and m 1 , respectively, and infected ticks with the invasive strain and co-infected (i.e., infected with both strains V 1 and V 2 ) ticks, m 2 and m c , respectively.Figure 3a shows the mathematical model and the routes of viral transmission considered between  6).Model parameters have been chosen as discussed in Section 2.2 and units for α 1 , β 1 and γ 1 as in Table 1.On the left, α 1 = 10 −6 , in the middle α 1 = 2 × 10 −5 , and on the right α 1 = 10 −4 .The parameters γ 1 and β 1 are varied along the x-axis and y-axis, respectively, from 0 to their maximum value listed in Table 1.Black curves represent the critical value R 0 = 1.
different tick compartments.The model corresponds to the following system of ODEs: where we have introduced We have assumed that the m 1 (m 2 ) population has transmission parameter α 1 (α 2 ) for V 1 (V 2 ), and the m c (co-infected) population has transmission parameter δ 1 for V 1 and δ 2 for V 2 , respectively.We have also slightly abused notation by writing Φ T = Φ.(8).Transmission rates are defined in Eq. ( 9).(b) Withinhost mathematical model of co-feeding transmission defined by Eq. (12).Transmission rates are defined in Eq. ( 9).(c) Alizon's (generalised) proposal for co-infection and cotransmission described in Eq. (13).Transmission rates are defined in Eq. ( 14).(d) Two-slot mathematical model of co-infection and co-transmission defined in Eq. (15).Transmission rates are defined in Eq. ( 14) and Eq. ( 16).

Basic reproduction number
The mathematical model defined by the system of ODEs ( 8) has a virusfree equilibrium (VFE), (m ⋆ 0 , 0, 0, 0), with m ⋆ 0 = Φ ν 0 .To compute its basic reproduction number, we make use of the next-generation matrix method, as illustrated in detail in Section 2.1.The T and V matrices are given by Thus, by computing the eigenvalues of the next-generation matrix, K = T (−V ) −1 , the basic reproduction number of system ( 8) can be shown to be Following the results from Ref. [18, Proposition 2.1], we can explore the boundary equilibria of system (8): 1.The virus-free equilibrium, E 0 = (m ⋆ 0 , 0, 0, 0), always exists.

Invasion reproduction number
We now assume that R 1 > 1, so that the endemic equilibrium E 1 of (8) exists.We write with E 1 = (m 0 , m 1 , 0, 0).We aim to calculate the invasion reproduction number of V 2 by means of the next-generation matrix method.To this end, we identify the invasive sub-system of V 2 of Eq. ( 8), linearise it around E 1 , compute its Jacobian matrix, and define the 2 × 2 matrices T and V .We can write , and The next-generation matrix, K = T (−V ) −1 , is given by with the type reproduction numbers R 22 , R 2c , R c2 , and R cc given by The eigenvalues of K are solutions of the following quadratic equation The invasion reproduction number, R I , is the largest eigenvalue of K, i.e., able to invade a tick population where the resident strain V 1 is endemic.

Alternative neutral models of co-feeding, co-infection, and cotransmission
The invasion reproduction number of the mathematical model from Section 3.2 is not neutral [23,37].By neutrality, we mean the following: in the limit when the invasive strain tends to the resident one, there should be no advantage for either strain, and thus, R I → 1.One can show for R I given by Eq. ( 10) that R I ̸→ 1.In fact, we have R 1 , for ϵ c = 0 and ϵ c = 1, respectively, under the assumption that infection does not affect the death rate, i.e., ν 0 = ν 1 = ν 2 = ν c .The issue of neutrality in co-infection models was brought up by Samuel Alizon in Ref. [37] and Lipsitch et al. in Ref. [35].We now present five alternative neutral formulations of the previous model.The first (and less optimal) option for obtaining a neutral model is to force R I → 1 and in turn, consider the constraints this condition imposes on some of the model parameters.The second one, as proposed by us to Samuel Alizon in private communication, is to consider a normalised invasion reproduction number; that is, define where by lim 2→1 R I , we mean the value of the invasion reproduction number in the limit when the invasive strain tends to the resident one (see Section 4.1).The third one generalises the mathematical model ( 8) by introducing the idea of within-host probability of invasion (see Section 4.2).The fourth one, as proposed by Alizon [37], is to consider a more general class of models, with doubly infected individuals (see Section 4.3).A final one that we propose in Section 4.4, is a generalisation of the approach of Alizon [23,37], which clearly articulates the issue of co-transmission.

A normalised invasion reproduction number
The invasion reproduction number given by Eq. ( 10) is not neutral.Let us then define a normalised invasion reproduction number, R N I , as follows where lim 2→1 means ν 2 → ν 1 , α 2 → α 1 , and δ 2 → δ 1 for our co-infection model (see Eq. ( 8)).So defined, it is clear that lim 2→1 R N I = 1, which is the desired neutrality condition.We note that the condition for the invasive strain to have the potential to become established is R I > 1.Now that we have introduced a normalised invasion reproduction number, this condition becomes

A model with within-host invasion
The fitness advantage of the invasive strain in model (8) stems from the assumption that V 2 can infect susceptible ticks and infected ticks by the resident strain with the same rate.However, this may not be realistic.For instance, a small amount of transmitted (invasive) virus may be less likely to establish infection in a tick that already has a high resident viral load, compared to a fully susceptible tick.The probability of within-host invasion will, thus, depend on the relative within-host fitnesses of the invasive and resident strains.Therefore, we can adapt the previous model (Eq.( 8)) by introducing the parameter ϕ i|j , which is the probability that strain i can establish co-infection in a tick already infected by strain j, given that there is transmission of strain i via co-feeding.This is similar to the super-infection framework described by Alizon in Ref. [36].The model can then be described by the following system of ODEs: where λ 1 , λ 2 and λ c are defined by Eq. ( 9).The transmission events of this model are summarised in Fig. 3b.In Appendix B we show that the invasion reproduction number for this model satisfies the desired neutrality condition.

A generalisation of Alizon's proposal
The mathematical model proposed by Alizon in Ref. [37] to obtain a neutral invasion reproduction number requires two additional populations (see Fig. 3c), namely the populations of doubly infected ticks with either V 1 or V 2 , denoted by M 1 and M 2 , respectively.Thus, there are six different tick compartments: m 0 , susceptible ticks, m 1 , m 2 , ticks (singly) infected with either V 1 or V 2 , M 1 , M 2 , doubly infected ticks with either V 1 or V 2 , and m c , co-infected ticks with both V 1 and V 2 .We note that the co-infection models in Ref. [23] do not consider co-transmission, but it is discussed in Ref. [37].Thus, in what follows, and to develop a mathematical model of co-infection and co-transmission in co-feeding ticks, we explain in detail what happens when a co-infected tick transmits virus to a singly infected tick.If cotransmission of both, resident and invasive, strains occurs, the singly infected tick can only acquire one new viral strain, since the mathematical model does not accommodate triply infected ticks.Therefore, if co-transmission takes place, a singly infected tick will acquire V 1 with probability δ 1 δ 1 +δ 2 , or V 2 with probability δ 2 δ 1 +δ 2 .Hence, the overall rate at which a co-infected tick transmits V 1 to a singly infected tick is (1 where the first term represents transmission of V 1 if no co-transmission, and the second term represents transmission of V 1 in the event of co-transmission.Similarly, the overall rate a co-infected tick transmits V 2 to a singly infected tick is δ 2 .We now write down the system of ODEs for Alizon's generalised mathematical model of a co-feeding tick population, with two circulating viral strains, which allows for co-infection and co-tranmission, and at most doubly infected ticks (with the same strain M 1 and M 2 , or with different ones m c ).We have where we define with ϵ c , the probability of co-transmission of the two viral strains by a coinfected (m c ) tick, and with ϵ 1 (ϵ 2 ) the probability of co-transmission by a doubly infected tick M 1 (M 2 ), respectively.The original model of co-infection with co-transmission defined in Ref. [37] assumed ϵ 1 = ϵ 2 = ϵ c = ϵ.We warn the reader that we have defined λ 1 and λ 2 to mean two different things in Eq. ( 14) and Eq. ( 9).We shall always clarify in what follows, which of the two definitions is implied.Fig. 3c shows the transmission events described by Eq. ( 13).We note that co-transmission by the M 1 tick population to a susceptible tick implies double transmission of the resident viral strain V 1 , and that co-transmission by the M 2 tick population (to a susceptible tick) implies double transmission of the resident strain V 2 .Finally, the parameter κ 1 (κ 2 ) is the rate of transmission of a single copy of V 1 (V 2 ) from an M 1 (M 2 ) tick to a susceptible one; thus, the factor of 2 in the previous expression for Λ 1 (Λ 2 ).We remind the reader that the model defined above includes death, immigration and transmission events.We have assumed each tick compartment has a different death rate, and immigration replenishes the susceptible tick compartment.In Appendix C we carefully derive the invasion reproduction number of this model and show its neutrality.

Two-slot model of co-feeding, co-infection and co-transmission
In the model defined by Eq. ( 13), a co-transmission event from a coinfected tick, in the m c compartment, to a susceptible tick implies the transmission of both viral strains at once.Here we extend the previous model to allow for the possibility that such a co-transmission event could instead result in the transmission of two copies of V 1 or two copies of V 2 .The idea of this generalised two-slot model is as follows: since ticks can be at most doubly infected, we assume each tick has two infection slots that can be occupied (or not).In the previous model, "co-transmission" to a susceptible tick meant transmission of both viral strains.In this model "co-transmission" means occupying both slots, in such a way, that the slots can be occupied by two copies of the same virus (leading to M 1 or M 2 ticks), or two different strains (leading to m c ticks).The dynamics of the two-slot model can be written as where λ 1 , λ 2 , λ 1,c , λ 1,c , Λ 1 , and Λ 2 have been defined in Eq. ( 13), and with Λ c given by In Appendix D we describe in great detail the transmission events considered in the two-slot mathematical model, show the existence of an endemic equilibrium for V 1 , and prove that the model leads to a neutral invasion reproduction number.

Numerical study of the invasion reproduction number
In Section 3 we have defined and computed the invasion reproduction number for a mathematical model of co-infection and co-transmission in cofeeding ticks.We have argued that such a model is not neutral, and have in turn proposed different mathematical models which do not suffer from such problem.We now propose a numerical study of the invasion reproduction number for the "not-neutral" model introduced in Section 3, as well as the invasion reproduction number for the model solutions proposed above to guarantee neutrality.
In what follows we assume that κ 1 = α 1 /2, and κ 2 = α 2 /2, which are appropriate choices when considering viral infections or micro-parasites [37].As discussed by Alizon in Ref. [37], if an infected host (by a certain viral strain) is re-infected by the exact same strain, we do not expect to see a change in its viral load (and hence in transmission rate).For co-infected ticks in the m c compartment (and infected by both V 1 and V 2 ), it is reasonable to hypothesise that potential within-tick interactions between the two strains do not lead to a change in transmission rates, when compared to doubly infected ticks in the M 1 or M 2 compartments.Thus, we assume δ 1 = κ 1 and δ 2 = κ 2 .Finally, and as justified earlier, we set ν 0 = ν 1 = ν 2 = ν c = υ 1 = υ 2 = 10 −2 per day.We also fix the immigration rate to be Φ = 2 ticks per day, and α 1 = 10 −4 per tick per day.These choices lead to a basic reproduction number of R 1 = 2 for the resident strain, V 1 .
In Fig. 4 we compare how different values of the transmission parameter, α 2 , and the co-transmission probability, ϵ c , affect the invasion reproduction number, R I , computed in the "not-neutral" scenario (10) (panel (a)), in the normalized proposal of Eq. ( 11) (panel (a ⋆ )), in the within-host model (12) (panel (b)), in Alizon's model with co-transmission (13) (panel (c)), and in the two-slot mathematical model (15) (panel (d)).In particular, ϵ c is varied along the x-axis, whereas the ratio α 2 /α 1 is varied from 0.5 to 1.5 along the yaxis.Black lines mark the contours where the invasion reproduction number, R I , is equal to 1.For Alizon's model and the two-slot extension, we have set ϵ 1 = ϵ 2 = 1 /2.For the within-host model, we define the probability that strain i can establish co-infection in a tick already infected by strain j as follows We note that the highest values of the invasion reproduction number occur in panel (a) and (b) of Fig. 4. In panel (a), the invasion reproduction number is clearly not neutral, since R I > 1 when α 2 = α 1 , and also for some regions of parameter space with α 2 < α 1 .For this model, if infected ticks with the invasive strain, V 2 , are rare compared to ticks infected with the resident strain, V 1 , then V 2 has an initial advantage over V 1 .Each tick infected with the invasive strain has the opportunity to infect a much larger number of ticks (m 0 + m 1 ), than those which can be infected by a tick from the m 1 compartment.This allows V 2 to invade the V 1 endemic system, for large enough values of α 2 and ϵ c .The co-transmission probability, ϵ c , affects the value of R I , since it changes the rate at which co-infected ticks transmit V 1 and V 2 to susceptible ticks, m 0 .These rates are δ 1 + ϵ c δ 2 and δ 2 + ϵ c δ 1 , respectively.Therefore a higher probability of co-transmission enables both strains to be transmitted more often.For the normalised invasion reproduction number in panel (a ⋆ ), R I = 1 when α 2 = α 1 , for every value of ϵ c , given its definition.When α 2 ̸ = α 1 , R I does depend on ϵ c , but less so than for the model of panel (a), since lim 2→1 R I increases with ϵ c .In panel (c), showing the invasion reproduction number for Alizon's model, when ϵ c = 1 /2 (equal to ϵ 1 and ϵ 2 ), R I = 1 for α 2 = α 1 .As ϵ c increases, so does the value of R I , since a higher co-transmission probability enables V 2 to be transmitted along with V 1 more often.The invasion reproduction number of the two-slot model in panel (d) behaves in a similar fashion.However, increasing ϵ c does not give as much of an advantage to the invasive strain, since co-transmission events can result in the transmission of two copies of V 1 .

Discussion and conclusions
In this paper, we consider the role of different transmission routes for a single vector-borne virus in a population of ticks and vertebrate hosts.We then study co-infection and co-transmission of two circulating vector-borne viral strains in a population of co-feeding ticks.We define and compute both the basic reproduction number and the invasion reproduction number, which provides the conditions under which a new variant can emerge (possibly endogenously from genomic reassortment).We illustrate how a classic and intuitive model of invasion was not, in fact, neutral with respect to the invading strain; that is, using this model to understand, for example, the minimum selective advantage that needs to be present for a invading strain to take hold of an endemic population (with the resident one) will privilege one strain over the other.This is not a problem per se, as it might be the correct model from a mechanistic perspective.However, it is important to characterise the underlying properties of a mathematical model, especially if it is intended to be used as part of an inference procedure.We also presented several alternative formulations of co-infection and co-transmission models that are, by definition neutral.We have shown that each model has distinct and specific behaviour concerning the invasion reproduction number.The take-home message of this review is that the assumptions used to model these important and complex infection systems matter, specially when making inferences about pathogens of potential pandemic emergence.In the real world, the choice of model, from the different alternatives presented and discussed here, will clearly depend on the virus, as well as the immunology and ecology of the hosts those viruses infect.
In conclusion, we note that while we have focused on deterministic models of tick-borne disease transmission, stochastic analogous may be considered instead, particularly when studying the invasion potential of a rare circulating viral strain [15,51,52].In a stochastic framework, the reproduction number is defined as a random variable rather than as an average [53], since its distribution encodes the probability of an epidemic occurring if a pathogen is introduced into a fully susceptible population by a small number of infected individuals [53].Thus, future work should include a study of the invasion reproduction number probability distribution, as well as an exploration of the issue of non-neutrality making use of stochastic approaches [54].Finally, given recent reports which indicate an increase in the number of Zika and Dengue virus co-infection cases in expanding co-endemic regions [55], it is of utmost importance to have suitable within-host mathematical models to study the impact of co-infection on viral infection dynamics. where When considering neutrality (i.e., the invasive strain is the same as the resident strain), all infected tick populations (m 1 , m 2 , and m c ) are infected with the same viral strain.Thus, we have ν 2 = ν c = ν 1 and α 2 = α 1 .We also set δ 1 = δ 2 = α 1 2 , representing that ticks in the m c compartment will transmit virus at the same overall rate as ticks in the m 1 compartment (i.e., δ 1 + δ 2 = α 1 ).Furthermore, we would expect ϕ 1|2 = ϕ 2|1 = 0, since in the within-host environment (a tick) the transmitted strain is likely to be rare compared to the established strain and if the resident strain is the same as the invasive strain, then both strains have the same within-host fitness, implying that the rare transmitted strain will have no within-host advantage over the established strain, and will be unable to establish co-infection in the host (tick).With these limits, the elements of the next-generation matrix become The eigenvalues are then b 11 and b 22 , which are equal to 1 and ϵ c , respectively, since m 0 = ν 1 α 1 .Thus, we can conclude that R I = 1, given ϵ c ∈ [0, 1].

Appendix C. Alizon's proposal
Alizon [37] proposed a model with doubly infected hosts (with the same viral strain), and which seemed a sufficient approach to achieve neutral-ity.We have considered a generalisation of the model originally proposed in Ref. [37], and which is described by the ODE system (13).By setting = 0, and m 2 = m c = M 2 = 0, one obtains the endemic equilibrium for the resident strain E 1 = ( m 0 , m 1 , 0, 0, M 1 , 0).We then compute the invasion reproduction number of the invasive strain by considering the invasive sub-system, linearised around the resident strain endemic equilibrium, E 1 : The Jacobian matrix of the invasive sub-system is and can be decomposed as follows Finally, the next-generation matrix, hosts.The assumption that κ = α/2 is realistic because we are considering viral infection.As Alizon mentioned in Ref. [37], if a host infected by a given strain is re-infected by the exact same strain, we do not expect to see a change in viral load (and hence in transmission rate); that is, a singly infected tick can become doubly infected with the same strain in this model, but becoming doubly infected does not affect its viral load.Under these conditions, the endemic equilibrium for the resident strain satisfies with Hence, we have • Transmission from a singly infected tick to a singly infected tick: • Transmission from a doubly infected tick to a susceptible tick: • Transmission from a co-infected tick to a susceptible tick: • Transmission from a co-infected tick to a singly infected tick: • Transmission from a doubly infected tick to a singly infected tick:

Appendix D.3. Invasion reproduction number
We can now compute the invasion reproduction number, R ′ I , of V 2 for the invasion sub-system, linearised around the endemic equilibrium E ′ 1 .We have the following Jacobian matrix: , which can be decomposed as follows We can compute the next-generation matrix, K ′ = T ′ [−V ′ ] −1 , given by: where We can obtain the invasion reproduction number, R ′ I , as a function of m ′ 0 by substituting Eq. (D.6) into K ′ and computing its largest eigenvalue.
carried out by means of a branching process approximation to calculate the probability that the new strain becomes established.5. Allen et al. [24] formulated a general epidemiological model for one vector species and one plant species with potential co-infection in the host plant by two viruses.First, the basic reproduction number is derived, and thus, conditions for successful invasion of a single virus are determined.Then, a new invasion threshold is derived to provide conditions for successful invasion of a second virus.6. White et al. [25] proposed a mathematical model for a two-pathogen, one-tick, one-host system with the aim of determining how long an invading pathogen persists within a tick population in which a resident pathogen is already established.7. Rovenolt et al. [29] developed a model of two co-infected host species to understand under which conditions co-infection can interfere with parasite-mediated apparent competition among hosts.8. Le et al. [30] studied a Susceptible-Infected-Susceptible (SIS) compartmental model with two strains and co-infection/co-colonization, incorporating five strain fitness dimensions under the same framework to understand coexistence and competition mechanisms.9. Alizon [36] discussed how multiple infections have been modelled in evolutionary epidemiology, presenting within-host models, super-infection frameworks, co-infection models, and some perspectives for the study of multiple infections in evolutionary epidemiology.In particular, he showed that a widely used co-infection model is not neutral as it confers a frequency-dependent advantage to rare neutral mutants.10.Alizon [37] studied the effect of co-transmission on virulence evolution when parasites compete for host resources.11.Bhowmick et al. [38] developed a compartment-based non-linear ordinary differential equation system to model the disease transmission cycle including blood-sucking ticks, livestock and humans.Sensitivity analysis of the basic reproduction number shows that decreasing the tick survival time is an efficient method to control the disease.They concluded that in the case of CCHFV transmission due to co-feeding, as well as trans-stadial and trans-ovarial transmission, are important routes to sustain the disease cycle.12. Hoch et al. [44] proposed a dynamic mechanistic model that takes into account the major processes involved in tick population dynamics and pathogen transmission with the aim of testing potential scenarios for pathogen control.13.Johnstone et al. [48] derived expressions for the basic reproduction number and the related tick type-reproduction number accounting for the observation that larval and nymphal ticks tend to aggregate on the same minority of hosts (tick co-aggregation and co-feeding).The pattern of tick blood meals is represented as a directed, acyclic, bipartite contact network.14.Belluccini [51] proposed both deterministic and stochastic models of coinfection with tick-borne viruses to investigate the role that different routes of transmission play in the spread of infectious diseases and to study the probability and timescale of co-infection events.15.Maliyoni et al. [52] investigated the impact of between-patch migration on the dynamics of a tick-borne disease on disease extinction and persistence making use of a system of stochastic differential equations.16.Lin et al. [55] studied the impact of Zika and Dengue virus co-infection on viral infection, examining viral replication activity in cells infected simultaneously, or sequentially.

1 Figure 1 :
Figure 1: A) If two viral strains, V A and V B , are co-circulating, the target cells, T , of an infected host will become singly infected (I A and I B ), and potentially co-infected (I C ).Co-infected cells have the potential to generate new viral progeny, different from that of the parental strains:V C ̸ = V A and V C ̸ = V B .B)A co-infected cell can lead to reassortment events, and produce up to 2 3 different reassortants.

Figure 3 :
Figure3: Illustrative diagrams of the mathematical models discussed in the paper for a population of co-feeding ticks with two viral strains.(a) Mathematical model of co-feeding transmission defined by Eq.(8).Transmission rates are defined in Eq. (9).(b) Withinhost mathematical model of co-feeding transmission defined by Eq.(12).Transmission rates are defined in Eq. (9).(c) Alizon's (generalised) proposal for co-infection and cotransmission described in Eq.(13).Transmission rates are defined in Eq. (14).(d) Two-slot mathematical model of co-infection and co-transmission defined in Eq.(15).Transmission rates are defined in Eq. (14) and Eq.(16).